### load data
#CHANGE working directory to wherever you put the data files
setwd("C:\\Users\\Jesse\\Desktop\\CNCopencode")
RE<- read.csv("individualData.csv")
TPss<- read.csv("trialData.csv")
SI<-read.csv("solutions1.csv")
SI2<-read.csv("solutions2.csv")


#define colors:
colPal<- c("#B4C2D5","#759EB8", "#729197", "#D15E5E","#663F46") 

suppressWarnings(suppressMessages(require(stargazer)))
suppressWarnings(suppressMessages(require(ggplot2)))
suppressWarnings( suppressMessages(require(lme4)))
suppressWarnings( suppressMessages(require(Matrix)))
library(dplyr)
library(lubridate)
library(quantreg)
library(survival)
library(survminer)
library(lme4)
library(nnet)
library(lqmm)


## define some convenience functions:
super.cluster.fun<-function(model, cluster)
{
  suppressWarnings(suppressMessages(require(multiwayvcov)))
  suppressWarnings(suppressMessages(require(lmtest)))
  vcovCL<-cluster.vcov(model, cluster)
  
  coef<-coeftest(model, vcovCL)
  w<-waldtest(model, vcov = vcovCL, test = "F")
  ci<-get_confint(model, vcovCL)
  
  return(list(coef, w, ci))
}
get_confint<-function(model, vcovCL){
  t<-qt(.975, model$df.residual)
  ct<-coeftest(model, vcovCL)
  est<-cbind(ct[,1], ct[,1]-t*ct[,2], ct[,1]+t*ct[,2])
  colnames(est)<-c("Estimate","LowerCI","UpperCI")
  return(est)
}
library(gridExtra)
library(grid)

grid_arrange_shared_legend <- function(..., nrow = 1, ncol = length(list(...)), position = c("bottom", "right")) {
  
  plots <- list(...)
  position <- match.arg(position)
  g <- ggplotGrob(plots[[1]] + theme(legend.position = position))$grobs
  legend <- g[[which(sapply(g, function(x) x$name) == "guide-box")]]
  lheight <- sum(legend$height)
  lwidth <- sum(legend$width)
  gl <- lapply(plots, function(x) x + theme(legend.position = "none"))
  gl <- c(gl, nrow = nrow, ncol = ncol)
  
  combined <- switch(position,
                     "bottom" = arrangeGrob(do.call(arrangeGrob, gl),
                                            legend,
                                            ncol = 1,
                                            heights = unit.c(unit(1, "npc") - lheight, lheight)),
                     "right" = arrangeGrob(do.call(arrangeGrob, gl),
                                           legend,
                                           ncol = 2,
                                           widths = unit.c(unit(1, "npc") - lwidth, lwidth)))
  grid.newpage()
  grid.draw(combined)
  
}


descTab<- array(0,c(4,7))

colnames(descTab)<- unique(TPss$network)
rownames(descTab) <- c("N", "mean(# correct)", "sd(# correct)","median(# correct)")

for(net in colnames(descTab)){
  descTab[1,net] <- length( TPss$numCorrectAtFinish[TPss$network ==net])
  descTab[2,net] <- mean( TPss$numCorrectAtFinish[ which(TPss$network ==net)])
  descTab[4,net] <- median( TPss$numCorrectAtFinish[ which(TPss$network ==net)])
  descTab[3,net] <- sd( TPss$numCorrectAtFinish[ which(TPss$network ==net)])
  
}
descTab
descTabLatex <- stargazer(descTab, title="Descriptive statistics by network treatment", label = "Tab:descriptives", font.size ="small" )


######################
#main results tables #
######################
#estimate models in main results:
library(lqmm) #for quantile regression models
TPss$network<- relevel(TPss$network, ref="HubAndSpoke")
q25HAS<- (lqm.counts(numCorrectAtFinish~network, data=TPss, M=200, tau=.25)) #for 25th quantile, relative to hub and spoke
q50HAS <- (lqm.counts(numCorrectAtFinish~network, data=TPss, M=200, tau=.5, control=lqmControl(loop_max_iter = 5000,loop_tol_ll = 2e-5)))
q75HAS <- (lqm.counts(numCorrectAtFinish~network, data=TPss, M=200, tau=.75))
mHAS  <- (glm(numCorrectAtFinish~network, data=TPss, family=quasipoisson)) #for mean, relative to hub and spoke

TPss$network<- relevel(TPss$network, ref="CorePeriphery")
q25CP<- (lqm.counts(numCorrectAtFinish~network, data=TPss, M=200, tau=.25))
q50CP <- (lqm.counts(numCorrectAtFinish~network, data=TPss, M=200, tau=.5, control=lqmControl(loop_max_iter = 5000,loop_tol_ll = 2e-5)))
q75CP <- (lqm.counts(numCorrectAtFinish~network, data=TPss, M=200, tau=.75))

mCP  <- (glm(numCorrectAtFinish~network, data=TPss, family=quasipoisson))

#latex output:
stargazer(mHAS,mHAS,mHAS, mHAS,mCP,mCP,mCP,mCP,
          
          se=list(q25HAS$tTable[,2], q50HAS$tTable[,2], NULL, q75HAS$tTable[,2],  q25CP$tTable[,2], q50CP$tTable[,2],NULL,  q75CP$tTable[,2]),
          coef=list(q25HAS$tTable[,1], q50HAS$tTable[,1], NULL, q75HAS$tTable[,1],  q25CP$tTable[,1], q50CP$tTable[,1],NULL, q75CP$tTable[,1]),
          title = "Number of people correct at end of trial",
          omit.stat=c("ser","f", "ll", "aic", "bic"), font.size="small", no.space=TRUE, model.names = FALSE, 
          dep.var.labels = rep("",8),
          label="tab:numCorrectQRresults",
          column.labels = c("25$^{th}$ \\%$^{ile}$", "median", "mean", "75$^{th}$ \\%$^{ile}$","25$^{th}$ \\%$^{ile}$", "median","mean", "75$^{th}$ \\%$^{ile}$" ))#,
          

######################################
#########visualization of CDFs:#######
######################################

legxloc<-.75 #for convenience, define a location to put legends
grid.arrange(
  
  ggplot(TPss[TPss$network %in% c( "CompleteClique", "SyntheticIsolates"),], aes( numCorrectAtFinish/(9), color=factor(network, levels=c("SyntheticIsolates", "CompleteClique", "LocallyClustered", "HubAndSpoke", "CorePeriphery", "OutCorePeriphery", "InCorePeriphery")) ) )+stat_ecdf(size=2)+theme_bw()+
    annotate("segment",xend = .135, yend = .35, x = .3, y = .35, colour = "black",arrow = arrow(length = unit(0.3,"cm")))+
    annotate("segment",x = .48, yend = .72, xend = .73, y = .72, colour = "black",arrow = arrow(length = unit(0.3,"cm")))+
    annotate("text",x = .22, y = .4, colour = "black", label="herding")+
    annotate("text",x = .6,  y = .77, colour = "black", label="learning")+ 
    theme(legend.position = c(legxloc, 0.2))+
    xlab("number correct")+
    ylab("cumulative percent of trials")+ 
    theme(legend.title = element_blank())+
    scale_color_manual(values=colPal[1:2],labels = c("Isolates", "Complete clique"))+ggtitle("A"),
  
  ggplot(TPss[TPss$network %in% c( "CompleteClique", "SyntheticIsolates","LocallyClustered"),], aes(numCorrectAtFinish/(9), color=factor(network, levels=c("SyntheticIsolates", "CompleteClique", "LocallyClustered", "HubAndSpoke", "CorePeriphery", "OutCorePeriphery", "InCorePeriphery"))))+stat_ecdf(size=2)+theme_bw()+
    theme(legend.position = c(legxloc, 0.2))+
    xlab("number correct")+
    ylab("cumulative percent of trials")+ 
    theme(legend.title = element_blank())+
    scale_color_manual(values=colPal[c(1,2,3)],labels = c("Isolates", "Complete clique", "Locally-Clustered"))+ggtitle("B"),
  
  ggplot(TPss[TPss$network %in% c("HubAndSpoke", "CompleteClique", "SyntheticIsolates"),], 
         aes(numCorrectAtFinish/(9), 
             color=factor(network, levels=c("SyntheticIsolates", "CompleteClique", "LocallyClustered", "HubAndSpoke", "CorePeriphery", "OutCorePeriphery", "InCorePeriphery"))))+
    stat_ecdf(size=2)+theme_bw()+
    theme(legend.position = c(legxloc, 0.2))+xlab("number correct")+
    ylab("cumulative percent of trials")+ 
    theme(legend.title = element_blank())+
    scale_color_manual(values=colPal[c(1,2,4)],labels = c("Isolates", "Complete clique", "Hub-and-spoke"))+ggtitle("C"),
  
  ggplot(TPss[TPss$network %in% c( "CompleteClique", "SyntheticIsolates","CorePeriphery"),], 
         aes(numCorrectAtFinish/(9), 
                color=factor(network, levels=c("SyntheticIsolates", "CompleteClique", "LocallyClustered", "HubAndSpoke", "CorePeriphery", "OutCorePeriphery", "InCorePeriphery"))))+stat_ecdf(size=2)+theme_bw()+
      theme(legend.position = c(legxloc, 0.2))+
    xlab("number correct")+
    ylab("cumulative percent of trials")+ 
    theme(legend.title = element_blank())+
    scale_color_manual(values=colPal[c(1,2,4)],labels = c("Isolates", "Complete clique", "Core-periphery"))+ggtitle("D"),
  
  ggplot(TPss[TPss$network %in% c( "CompleteClique", "SyntheticIsolates","OutCorePeriphery"),], aes(numCorrectAtFinish/(9), color=factor(network, levels=c("SyntheticIsolates", "CompleteClique", "LocallyClustered", "HubAndSpoke", "CorePeriphery", "OutCorePeriphery", "InCorePeriphery"))))+stat_ecdf(size=2)+theme_bw()+
      theme(legend.position = c(legxloc, 0.2))+xlab("number correct")+ylab("cumulative percent of trials")+ theme(legend.title = element_blank())+scale_color_manual(values=colPal[c(1,2,5)],labels = c("Isolates", "Complete clique", "Out Core-periphery"))+ggtitle("E"),
  
  ggplot(TPss[TPss$network %in% c("CompleteClique", "SyntheticIsolates","InCorePeriphery"),], 
         aes(numCorrectAtFinish/(9), color=factor(network, levels=c("SyntheticIsolates", "CompleteClique", "LocallyClustered", "HubAndSpoke", "CorePeriphery", "OutCorePeriphery", "InCorePeriphery"))))+stat_ecdf(size=2)+theme_bw()+
      theme(legend.position = c(legxloc, 0.2))+
    xlab("number correct")+ylab("cumulative percent of trials")+ theme(legend.title = element_blank())+scale_color_manual(values=colPal[c(1,2,5)],labels = c("Isolates", "Complete clique", "In core-periphery"))+ggtitle("F"),
  ncol=2
)



################################################################
# prob. of having same answers as neighbors in phase 3 of expt #
################################################################

SI$position <- ""
SI$position[which(SI$network == "OutCorePeriphery" & SI$node %in% 1:3)] <- "COcore"
SI$position[ which(SI$network == "OutCorePeriphery" & SI$node %in% c(4:9))]  <- "COperi"
SI$position[ which(SI$network == "InCorePeriphery" & SI$node %in% 1:3)] <- "CIcore"
SI$position[which(SI$network == "InCorePeriphery" & SI$node %in% c(4,6,8))] <- "CIperi1"
SI$position[which(SI$network == "InCorePeriphery" & SI$node %in% c(5,7,9))] <- "CIperi0"
SI$position[ which(SI$network == "CorePeriphery" & SI$node %in% c(1:3))] <- "CBcore"
SI$position[ which(SI$network == "CorePeriphery" & SI$node %in% c(4:9))] <- "CBperi"
SI$position[ which(SI$network == "LocallyClustered" )] <- "SW"
SI$position[ which(SI$network == "CompleteClique" )]  <- "nineClique"
SI$position[ which(SI$network == "HubAndSpoke" & SI$node == 1)]<- "C1core"
SI$position[ which(SI$network == "HubAndSpoke" & SI$node != 1)]<- "C1peri"

SI2$position <- ""
SI2$position[which(SI2$network == "OutCorePeriphery" & SI2$node %in% 1:3)] <- "COcore"
SI2$position[ which(SI2$network == "OutCorePeriphery" & SI2$node %in% c(4:9))]  <- "COperi"
SI2$position[ which(SI2$network == "InCorePeriphery" & SI2$node %in% 1:3)] <- "CIcore"
SI2$position[which(SI2$network == "InCorePeriphery" & SI2$node %in% c(4,6,8))] <- "CIperi1"
SI2$position[which(SI2$network == "InCorePeriphery" & SI2$node %in% c(5,7,9))] <- "CIperi0"
SI2$position[ which(SI2$network == "CorePeriphery" & SI2$node %in% c(1:3))] <- "CBcore"
SI2$position[ which(SI2$network == "CorePeriphery" & SI2$node %in% c(4:9))] <- "CBperi"
SI2$position[ which(SI2$network == "LocallyClustered" )] <- "SW"
SI2$position[ which(SI2$network == "CompleteClique" )]  <- "nineClique"
SI2$position[ which(SI2$network == "HubAndSpoke" & SI2$node == 1)]<- "C1core"
SI2$position[ which(SI2$network == "HubAndSpoke" & SI2$node != 1)]<- "C1peri"


SI2$solution<- as.factor(SI2$solution)

library(nnet)
SI$middlephase <- ifelse(SI$seconds>300 & SI$seconds<=600, 1, 0)
SI$lastphase <- ifelse(SI$seconds>600, 1, 0)
SI$shareCorrect <- ifelse(SI$altersWAnswer>0,SI$correctAlters/SI$altersWAnswer,0)
SI$shareRH <- ifelse(SI$altersWAnswer>0,SI$RHalters/SI$altersWAnswer,0)
SI$solution<- as.factor(SI$solution)
SI$solution<-relevel(SI$solution, ref="OA") 
SI$phase <- ifelse(SI$seconds<=300,1,ifelse(SI$seconds<=600,2,3)) 
SI$phase <- as.factor(SI$phase)

SI2$middlephase <- ifelse(SI2$seconds>300 & SI2$seconds<=600, 1, 0)
SI2$lastphase <- ifelse(SI2$seconds>600, 1, 0)
SI2$shareCorrect <- ifelse(SI2$altersWAnswer>0,SI2$correctAlters/SI2$altersWAnswer,0)
SI2$shareRH <- ifelse(SI2$altersWAnswer>0,SI2$RHalters/SI2$altersWAnswer,0)
SI2$solution<- as.factor(SI2$solution)
SI2$solution<-relevel(SI2$solution, ref="OA") 
SI2$phase <- ifelse(SI2$seconds<=300,1,ifelse(SI2$seconds<=600,2,3)) 
SI2$phase <- as.factor(SI2$phase)


mnlPos <- multinom(solution~(as.factor(correctAlters))*position, 
                   data=SI2[which(SI2$phase==3),],maxit = 1000)

newdataPos<- unique(SI2[which(SI2$phase==3),c("position", "correctAlters","network")])
newdataPos <- newdataPos[!newdataPos$network %in% c("", "SyntheticIsolates"),] 
newdataPos <- newdataPos[!newdataPos$position %in% c("CIperi0"),] 

newdataPos<-newdataPos[order(newdataPos$correctAlters),] 
newdataPos$CApreds<- predict(mnlPos, newdata=newdataPos, type='probs')[,2]
newdataPos$periphery <- grepl(newdataPos$position, pattern = "peri")

mnlPos2 <- multinom(solution~(as.factor(RHalters))*position, 
                    data=SI2[which(SI2$phase==3),],maxit = 1000)

newdataPos2 <- unique(SI2[which(SI2$phase==3),c("position", "RHalters","network")])
newdataPos2 <- newdataPos2[!newdataPos2$network %in% c("", "SyntheticIsolates"),] 
newdataPos2 <- newdataPos2[!newdataPos2$position %in% c("CIperi0"),] 

newdataPos2<-newdataPos2[order(newdataPos2$RHalters),] 
newdataPos2$RHpreds<- predict(mnlPos2, newdata=newdataPos2, type='probs')[,4]
newdataPos2$periphery <- grepl(newdataPos2$position, pattern = "peri")

newdataPos$network<- as.factor(newdataPos$network) 
newdataPos2$network<- as.factor(newdataPos2$network) 


grid.arrange(
  ggplot(newdataPos, aes(as.numeric(correctAlters), CApreds, linetype=periphery))+geom_path(size=1)+ xlab("number of correct visible answers")+ylab("probability of having correct answer registered in phase 3")+
    scale_x_continuous(limits=c(0,8), breaks = 0:8, minor_breaks = 0:8)+facet_wrap(~network)+
    scale_y_continuous(limits=c(0,1), breaks = seq(0,1,.1))+facet_wrap(~network)+
    theme_bw() +ggtitle("Probability of having answers in common with neighbors: correct answers")+theme(legend.position = "bottom"),
  
  
  ggplot(newdataPos2, aes(as.numeric(RHalters), RHpreds, linetype=periphery))+geom_path(size=1)+ xlab("number of red herring visible answers")+ylab("probability of having red herring registered in phase 3")+
    scale_x_continuous(limits=c(0,8), breaks = 0:8, minor_breaks = 0:8)+facet_wrap(~network)+
    scale_y_continuous(limits=c(0,1), breaks = seq(0,1,.1))+facet_wrap(~network)+
    theme_bw()+ggtitle("Probability of having answers in common with neighbors: red herrings")+theme(legend.position = "bottom"),
  ncol=1
) 


########################
 # FIGURE:

SI2$propRH<- ifelse(SI2$altersWAnswer==0,0,SI2$RHalters /SI2$altersWAnswer)
SI2$propCA<- ifelse(SI2$altersWAnswer==0,0,SI2$correctAlters /SI2$altersWAnswer)
SI2$propMatch<- ifelse(SI2$totalAlters==0,0,SI2$matchingAlters /SI2$totalAlters)
SI2$periphery<- grepl(SI2$position, pattern = "peri")


grid.arrange(
  ggplot(SI2[which(!SI2$network %in% c("", "synthNT")),], 
         aes(seconds, unqSolsVisible, linetype=periphery, color=periphery))+
    geom_smooth(se=FALSE,method = "loess", span=.15, method.args=list(degree=0) )+
    scale_x_continuous(limits=c(0,900), breaks = c(0,300,600,900))+facet_wrap(~network)+ggtitle("Unique solutions visible")+theme_bw()+theme(legend.position = "bottom")+scale_color_manual(values=colPal[c(2,3)]),
  
   ggplot(SI2[which(SI2$seconds>=840 &SI2$solution!="correct" & !SI2$network %in% c("", "synthNT")),], 
         aes(periphery,correctAlters, linetype=periphery, fill=periphery))+facet_wrap(~network)+
    geom_boxplot()+theme_bw()+theme(legend.position = "bottom")+ggtitle("Opportunities to learn")+scale_fill_manual(values=colPal[c(1,2)]) #geom_histogram(bins=9)
  
)
##############################

# MECHANISM C:
# descriptives for the order of events in terms of the core adopting the correct answer:

RE$position <- ""
RE$position[which(RE$network == "OutCorePeriphery" & RE$node %in% 1:3)] <- "COcore"
RE$position[ which(RE$network == "OutCorePeriphery" & RE$node %in% c(4:9))]  <- "COperi"
RE$position[ which(RE$network == "InCorePeriphery" & RE$node %in% 1:3)] <- "CIcore"
RE$position[which(RE$network == "InCorePeriphery" & RE$node %in% c(4,6,8))] <- "CIperi1"
RE$position[which(RE$network == "InCorePeriphery" & RE$node %in% c(5,7,9))] <- "CIperi0"
RE$position[ which(RE$network == "CorePeriphery" & RE$node %in% c(1:3))] <- "CBcore"
RE$position[ which(RE$network == "CorePeriphery" & RE$node %in% c(4:9))] <- "CBperi"
RE$position[ which(RE$network == "LocallyClustered" )] <- "SW"
RE$position[ which(RE$network == "CompleteClique" )]  <- "nineClique"
RE$position[ which(RE$network == "HubAndSpoke" & RE$node == 1)]<- "C1core"
RE$position[ which(RE$network == "HubAndSpoke" & RE$node != 1)]<- "C1peri"




RE$correctAtEndOrder<- 0 
RE$correctAtEndSeconds<- 0 

RE$periphery <- grepl(RE$position, pattern = "peri")


for(id in unique(RE$trial_id[which(RE$network !="SyntheticIsolates")])){
  if(any(RE$correctAtEnd[which(RE$trial_id ==id)]==1)){
    sirows<- SI[which(SI$trial_id==id),]
    sirows <- sirows %>% group_by(node) %>% filter(seconds == max(seconds) & solution == "correct")    
    for(nn in unique(sirows$node)){
      pos <- which(sirows$node == nn) 
      RE$correctAtEndOrder[ which(RE$trial_id == id & RE$node == nn)] <- pos
    }
  }
}

for(id in unique(RE$trial_id[which(RE$network !="SyntheticIsolates")])){
  if(any(RE$correctAtEnd[which(RE$trial_id ==id)]==1)){
    sirows<- SI[which(SI$trial_id==id),]
    sirows <- sirows %>% group_by(node) %>% filter(seconds == max(seconds) & solution == "correct")    
    for(nn in unique(sirows$node)){
      secs <- sirows$seconds[which(sirows$node == nn) ]
      RE$correctAtEndSeconds[ which(RE$trial_id == id & RE$node == nn)] <- secs
    }
  }
}

permarray<- array(0, c(10000,9))
trialstodo <- RE$trial_id[which(RE$network == "HubAndSpoke")]
for(i in 1:nrow(permarray)){
  tr<- sample(trialstodo, size=1) 
  permarray[i,]<- sample(RE$correctAtEndOrder[which(RE$trial_id == tr)])
  
}

coreOrderPerm<- array(0, c(10000, 10)) #holds the frequencies from permuted data for the core occupying each ordinal position 
colnames(coreOrderPerm)<- c(0:9)
for(i in 1:nrow(coreOrderPerm)){
  tab  <- table(sample(permarray[,1], 31))
  coreOrderPerm[i,names(tab)] <- tab
}

realCoreTab <- table(RE$correctAtEndOrder[which(RE$network == "HubAndSpoke" & RE$node==1)])



CPpermarray<- array(0, c(10000,9))
CPtrialstodo <- RE$trial_id[which(RE$network == "CorePeriphery")]
for(i in 1:nrow(CPpermarray)){
  tr<- sample(CPtrialstodo, size=1) 
  CPpermarray[i,]<- sample(RE$correctAtEndOrder[which(RE$trial_id == tr)])
  
}

CPcoreOrderPerm<- array(0, c(10000, 10)) #holds the frequencies from permuted data for the core occupying each ordinal position 
colnames(CPcoreOrderPerm)<- c(0:9)
for(i in 1:nrow(CPcoreOrderPerm)){
  theserows <- sample.int(nrow(CPcoreOrderPerm), 30)
  tab  <- table(CPpermarray[theserows,1:3])
  CPcoreOrderPerm[i,names(tab)] <- tab
}

CPrealCoreTab <- table(RE$correctAtEndOrder[which(RE$network == "CorePeriphery" & RE$node %in% 1:3)])

temp <- RE[which(RE$network == "HubAndSpoke" ),] 
temp <- temp[order(temp$correctAtEndSeconds),] %>%
  group_by(trial_id) %>% mutate(coreTime = correctAtEndSeconds[which(node == 1)])

HASGLM <- glm(correctAtEnd~I(coreTime == 0), data=temp[which(temp$node !=1 & (temp$correctAtEndSeconds ==0 | temp$coreTime< temp$correctAtEndSeconds)),], family="binomial")

temp <- RE[which(RE$network == "CorePeriphery" ),] 
temp <- temp[order(temp$correctAtEndSeconds),] %>%
  group_by(trial_id) %>% mutate(coreCorrect = sum(correctAtEndSeconds[which(node %in% 1:3)]!=0),
                                coreTime = max(correctAtEndSeconds[which(node %in% 1:3)])) 



CPGLM <- glm(correctAtEnd~as.factor(coreCorrect), 
             data=temp[which(temp$node %in% 4:9 & (temp$correctAtEndSeconds ==0 | temp$coreTime< temp$correctAtEndSeconds)),], family="binomial") 

# In the Following code, the \Sexpr{} expressions should be run to fill in the sentences with results (or compile with knitr)
##############################################################################################################################
# In the hub-and-spoke network, in cases in which the central node has already adopted the correct answer that they maintain through the end of the trial, 
# \Sexpr{100*round(predict(HASGLM, data.frame(coreTime = 1), type = "response"), digits=2)}\% of peripheral nodes that have not yet adopted the correct answer go on to do so, as compared to 
# \Sexpr{100*round(predict(HASGLM, data.frame(coreTime = 0), type = "response"), digits=2)}\% for cases in which the central node has not yet adopted the correct answer (the difference is statistically significant at $p= $
#                                                                                                                                                                           \Sexpr{round(summary(HASGLM)$coefficients[2,4], digits=6)}).  In Core-periphery networks, the probability that a peripheral node is correct at the end of the trial goes up sharply if two or three core nodes have already adopted the correct answer. The probability of completing the trial with the correct answer is 
# \Sexpr{100*round(predict(CPGLM, data.frame(coreCorrect = 0), type = "response"), digits=2)}\% if no core nodes have already adopted the correct answer, but 
# \Sexpr{100*round(predict(CPGLM, data.frame(coreCorrect = 2), type = "response"), digits=2)}\% if two core nodes have, and
# \Sexpr{100*round(predict(CPGLM, data.frame(coreCorrect = 3), type = "response"), digits=2)}\% if three core nodes have already adopted the correct answer ($p =$ 
#                                                                                                                                                              \Sexpr{round(summary(CPGLM)$coefficients[3,4], digits=3)} and 
#                                                                                                                                                            \Sexpr{formatC(format = "e",summary(CPGLM)$coefficients[4,4], digits=2)}, respectively). 

temp <- RE[which(RE$network == "HubAndSpoke"),] 
temp[order(temp$correctAtEndSeconds),] %>%
  group_by(trial_id) %>% mutate(coreTime = correctAtEndSeconds[which(node == 1)]) %>% 
  filter(correctAtEnd ==1) %>%
  ggplot(aes(correctAtEndSeconds, correctAtEndOrder, group=factor(trial_id)))+geom_step(alpha=.5)+
  geom_point(aes(color=node==1, shape=node==1), size=2)+
  facet_wrap(~ifelse(coreTime !=0, "core correct", "core incorrect"), ncol=1)+
  scale_y_continuous(breaks=1:9, minor_breaks = 1:9)+ylab("cumulative number correct")+
  scale_x_continuous(breaks=c(0,300,600,900))+xlab("time of final correct answer")+theme_bw()+theme(legend.position = "bottom")+
  scale_color_manual(name="Position: ", labels=c("periphery", "core"),values=colPal[c(2,4)])+guides(size = FALSE)+scale_shape_discrete(name="Position: ", labels=c("periphery", "core"))

